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I i ! . INTRODUCTION 

Projectiles  carrying  liquid  payloads  have  been  flown  for  at  least 
fifty  years1.  Quite  often,  their  flight  was  erratic  because  an  insta- 
bility developed  due  to  the  presence  of  the  liquid.  The  work  of 
Stewartson2  provided  an  understanding  of  this  phenomenon;  he  calculated 
the  natural  frequencies  of  .free  oscillation  of  a spinning  liquid  in  a 
cylindrical  container  and  demonstrated  that  a resonance  between  these 
frequencies  and  the  nutational  frequency  of  the  projectile  causes  the 
instability.  This  theory  and  important  modifications  to  it  by 
Wedemeyer3 4  have  provided  rational  design  methods  for  liquid-filled 
projectiles;  predictions 'ftom^these  have  been  verified  by  range  firings 
and  laboratory  experiments ; in;  which  the  conditions  for  the  theory  were 
satisfied.  • '••{..•  . y ' . . 

;v.d'  . • . = •* 

There  are  three  restrictive "assumptions  in  the  theory:  (i)  the 
amplitude  of  the  disturbance  to -the  .rotating  fluid  is  small;  (ii)  vis- 
cous effects  are  neglected;  (iii)  the  fluid  is  in  solid  body  rotation. 

An  adequate  treatment  of  large  amplitude  disturbances  is  not  yet  avail- 
able. Wedemeyer3  made  viscous veotfrect ions  to  Stewartson's  theory  which 
considerably  extended  its  usefulness^  Relaxing  restriction  (iii)  is 
the  objective  of  the  present-work..- 

The  time  required  for  a substantial  amount  of  fluid  to  achieve 
solid  body  rotation  is  called  the.ispin-mp  time,  t . If  t is  small 

compared  to  the  projectile  flight  time,  Stewartson's  assumption  (iii) 
is  reasonable,  and  his  theory  can  be  used.  The  applicability  of 
assumption  (iii)  depends  on  the  parameters  of  the  projectile-gun 
system.  When  (iii)  is  violated,  i't  is  necessary  to  calculate  the 
frequencies  of  the  liquid  durihg  ttie“'spin-up  phase  of  its  motion.  Pre- 
vious attempts  to  do  this  were  not  successful  because  of  the  complexity 
of  the  problem.  . . . , ... 

It  is  first  necessary  to  know ; the  basic  flow  during  spin-up.  The 
basic  theory  for  this  was  developed  by  Wedemeyer1*;  if  viscous  diffusion 

1 . Engineering  Design  ffgncffioQfc.,.  Liquid-Filled  Projectile  Design,  AMC 
Pamphlet  No.  70  6-165 , b S,^  Army'  Materiel  Development  and  Readiness 

Command,  Washington,  2?»  C.  Apvil  1969.  AD  855719. 

2.  K.  Stewartson,  " On  the  Stability  of  a Spinning  Top  Containing 
Liquid ,"  J . Fluid  Me ah. , Vol.  5,  No.  4,  September  1959,  pp.  577-592. 

5.  E.  H.  Wedemeyer,  "Visoous.  Corrections  to  Stewartson's  Stability 

Criterion,"  BRL  Report -No.  12  2S,  Aberdeen  Proving  Ground,  Maryland, 
June  1966.  AD  489687.  ' 

4.  E.  H.  Wedemeyer,  "The  Unsteady -.Flow  Within  a Spinning  Cylinder,"  J . 
Fluid  Mech.,  Vol.  20,.  Part  5,  19$ 4,  pp.  285-299.  Also,  see  BRL 
Report  No.  1252,  Aberdeen  Proving  Ground,  Maryland,  October  1962. 

AD  421846. 
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effects  are  neglected  in  his  model,  a simple  solution  for  the  spin-up 
motion  is  obtained.  Although  this  simple  solution  is  useful  for  many 
purposes,  it  is  inadequate  for  the  frequency  calculation.  We  have 
obtained  numerical  solutions  to  the  spin-up  problem,  including  viscous 
diffusion.  An  efficient  program  for  this  is  required  since  the  results 
furnish  input  to  the  frequency  computation.  The  latter  is  accomplished 
by  solving  an  eigenvalue  problem  for  a sixth-order  differential  system, 
the  same  system  of  perturbation  equations  one  would  obtain  in  performing 
a linear  stability  analysis  of  the  flow. 

The  Reynolds  numbers  of  interest  are  large,  or  Ekman  numbers 
are  small,  so  that  the  difficulties  encountered  in  solving  the 
eigenvalue  problem  for  the  Orr-Sommerfeld  equation  are  present  here 
also.  In  recent  years  several  techniques  have  been  developed  to  sur- 
mount these  difficulties;  we  choose  the  orthonormalization  technique  to 
solve  our  eigenvalue  problem.  Because  the  perturbation  equations  have 
a singularity  on  the  axis,  an  analytical  solution  must  be  obtained  there 
and  matched  with  the  numerical  solution. 

Although  the  problem  of  primary  interest  copcerns  spin-up.  Part  I 
of  this  report  contains  only  the  special  case  of  solid-body  rotation, 
because  the  authors  felt  it  important  to  give  a detailed  description  of 
the  numerical  procedures  yet  keep  the  report  a reasonable  length.  This 
division  is  feasible  because:  (i)  the  computational  method  is  general, 
applicable  both  to  solid-body  rotation  and  spin-up  flow;  (ii)  the 
presentation  will  not  be  hampered  by  complexities  which  arise  in  spin- 
UP  (e-g->  presence  of  Ekman  layers  and  critical  layers,  necessity  of 
obtaining  basic  flows  numerically) ; (iii)  checks  on  the  validity  of 
the  method  can  be  made  with  results  obtained  here.  Part  II  will  deal 
with  the  spin-up  period,  treated  in  a quasi-steady  manner  so  that  the 
same  computational  method  can  be  used. 

II.  GOVERNING  EQUATIONS  AND  BOUNDARY  CONDITIONS 

To  calculate  the  natural  frequencies  of  the  spinning  liquid  we 
assume  a small  disturbance  to  the  basic  flow.  The  Navier-Stokes 
equations  for  three-dimensional  incompressible  flow  are  linearized  to 
give  the  perturbation  equations.  If  the  basic  flow  is  solid-body 
rotation,  the  equations  can  be  found  in  Reference  3;  the  more  general 
form,  for  spin-up,  is  given  in  References  5 and  6.  The  non-dimensional 
flow  variables  are 


5.  Y.  M.  Lynn,  "Free  Oscillations  of  a Liquid  During  Spin-Up ,"  BRL 
Report  No.  1663,  Aberdeen  Proving  Ground,  Maryland,  August  1973. 

AD  769710 . 

6.  R.  Sedney  and  N.  Gerber,  "Perturbation  Equations  for  Liquid  Spin-Up 
in  Cylindrical  Cavities,"  BRL  Technical  Report  (in  preparation) . 
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(1) 


* * * _ * 
u = U + u,  v=V+v,  w=W+w,  p=P+p  , 

where  u,  v,  and  w are  total  radial,  azimuthal,  and  axial  components  of 
velocity,  respectively,  and  p is  the  pressure.*  U,  V,  W,  and  P are  the 
basic  unperturbed  variables,  and  u,  v,  w,  and  p are  the  small  perturba- 
tions in  the  sense  that  they  are  small  compared  to  U,  V,  W,  and  P. 

The  container  spin  rate  is  ft;  a is  the  container  radius  and  c is 
its  half-height;  v is  the  liquid  kinematic  viscosity.  Length,  velocity, 
pressure,  and  time  are  non-dimensionalized  by  a,  a ft,  pft2a2,  and  ft”  , 
respectively,  where  p is  the  liquid  density.  Dimensionless,  non- 
rotating, cylindrical  coordinates  r,  0,  and  z are  used,  with  z = 0 at 
one  endwall;  dimensionless  time  is  denoted  by  t.  For  solid-body  rota- 
tion (the  only  case  considered  here)  the  basic  flow  is  given  by 

U = 0,  V = r,  W = 0,  3P/3r  = r . (2) 

On  substituting  (1)  and  (2)  into  the  Navier-Stokes  equations  and 
linearizing,  we  obtain  the  perturbation  equations: 


* 

(ru)r 

* 

+ v0 

* 

+ r w 

z 

= 0 

(3a) 

* 

ut 

* 

+ U0 

* 

- 2v  = 

* 

- Pr  + Re 

1 (V2  u - 

r"2 

* 

u - 

2r"2  v0) 

(3b) 

* 

vt 

* 

+ v0 

* 

+ 2u  = 

-i  * 

- r P0  + 

Re"1  (V2 

* 

v - 

r"2 

* -r,  * 

v + 2 r 2 ufl) 

(3c, 

* 

wt 

* 

- w0 = 

* 

- p + Re 
r z 

1 V2  w 

(3d) 

where  subscripts  denote  partial  derivatives, 

V2  = 32/3r2  + r"1  3/3r  + r”2  3 2/30 2 + 32/3z2  , (4) 

and  the  Reynolds  number  is 


Re  = ft  a2/v  . (5) 

We  assume  that  the  disturbance  can  be  expressed  as  a superposition 
of  modes;  i.e.,  a triple  Fourier  expansion  of  the  disturbances  in  0,  z, 
and  t,  with  coefficients  functions  of  r: 
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where 


* 

u = 

A 

U 

(r) 

COS 

Kz 

exp 

[i 

(Ct  • 

- m0)] 

(6a) 

* 

V = 

V 

(r) 

cos 

Kz 

exp 

[i 

(Ct  - 

- me)] 

(6b) 

* 

w = 

w 

(r) 

sin 

Kz 

exp 

[i 

(Ct  ■ 

- me)] 

(6c) 

T3  si- 
ll 

A 

P 

(r) 

cos 

Kz 

exp 

[i 

(Ct  - 

- me)] 

(6d) 

K = 

kiT/  (2< 

z) 

(7) 

In  this  form  the  disturbances  are  complex,  the  real  parts  being  the 
physical  quantities;  also  the  functions  u(r),  v(r),  w(r),  and  p(r)  are 
complex.  The  integers  m = 0,  ±1,  and  k = 1,  2,  ...  are  azimuthal 

and  axial  wave  numbers,  respectively;  m expresses  the  condition  that  the 
disturbance  is  periodic  in  the  interval  0 < 0 < 2ir.  The  non-dimensional 
quantity  C = + i is  the  eigenvalue  we  seek  to  determine.  The 

natural  frequency  is  equal  to  C U and  the  decay  rate  of  the  disturbance 

is  equal  to  C^fi. 


* The  boundary  conditions  at  the  end  walls  (z  = 0,  2c)  are  u = v = 
w = 0.  The  form  assumed  in  (6)  allows  only  the  condition  of  zero 
velocity  normal  to  the  wall,  w = 0,  to  be  satisfied.  However,  the  modal 
decomposition,  or  separation  of  variables,  requires  the  trigonometric 
functions;  consequently,  the  no-slip  condition  on  the  endwalls,  u = v = 
0,  is  not  satisfied.  A boundary  layer  type  of  correction  is  made  to 
the  solution  to  account  for  this.  This  correction,  described  in 
Appendix  A,  is  analogous  to  that  made  by  Wedemeyer  at  both  the  sidewall 
and  endwalls.  The  correction  has  a significant  effect  on  . 

"k  "k  "k  k 

When  u,  v,  w,  and  p from  (6)  are  substituted  into  (3),  a set  of 
homogeneous  linear  ordinary  differential  equations  is  obtained  for  u, 
v,  vi,  p.  These  must  be  converted  to  canonical  form  to  be  integrated 
numerically;  i.e.. 


y!  = dyi/dr  = O,  Y1>  Y2>  •••  Yg) > 


i = 1 , 2 ...  6 , 


where 

yi  = u ^ ■ v’ 

A A A 

y2  = u - i v y4  = w 


y5  = w 


(8) 
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The  form  for  y is  convenient  for  satisfying  boundary  conditions.  After 
the  required  manipulations  are  performed  the  following  sixth  order 


system  of  equations  is  obtained: 

yj  = (m  - 1)  r-1  yx  - m r-1  y2  - K y4  (9a) 

y2  '=  (m  - 1)  r_1  y^^  - m r-1  y2  - i y3  - K y4  (9b) 

y^  = 2 (Re  + i mr-2)  y^  + i (B  + r"2) (y2  - y ^ - r-1  y3  - (9c) 

i m Re  r_1  y^ 

y'  •=  y (9d) 

y4  y5 

y3  = B y4  - r"1  y5  - K Re  y6  (9e) 


y^  = -Re” 1 B y4  + i Re-1  (2  Re  + i m r"2)  (y2  - + i ■ Re"1  r"1  y3  - 

(9f) 

K Re"1  y5  , 


where 

B = (m2/r2)  + K2  + i Re  (C  - m)  . 

The  terminal  conditions  at  the  sidewall  (r  = 1)  are  u (1)  = v (1)  = 
w (1)  = 0,  or 

y1  (1)  = y2  (1)  = y4  (1)  = 0 . (io) 

The  initial  conditions  depend  on  the  particular  problem  being  solved. 

For  a filled  cylinder  with  an  inner  concentric  rod,  or  central  burster, 
the  no- si ip  conditions  apply  at  the  inner  wall,  r = b: 

yl  (b-*  = y2  ^ = y4  = 0 • 
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For  a filled  cylinder  with  no  inner  rod  the  inner  boundary  is  r = 0, 
and  the  boundary  conditions6  depend  on  the  azimuthal  wave  number: 


m = 0: 

(0)  = y2 

(0)  = y5 

(0)  = 0 

(12a) 

m = 1 : 

y2 

t — \ 

O 
' — / 

II 

X 

4^ 

(°)  - y6 

ro)  = o 

(12b) 

m > 1 : 

yl 

(0)  = y2 

(0)  = y4 

«;o)  = o . 

(12c) 

Their  derivation  requires  only  kinematics  and  single-valuedness  (see 
Batchelor,  E.K.  and  Gill,  A.E.,  J.  Fluid  Mech,  14,  529). 

The  wave  numbers  m and  k are  assigned,  and  the  values  of  c and 
Re  are  given.  Then  (9),  (10),  and  (11)  or  (12)  constitute  an  eigen- 
value problem.  The  radial  modes  are  designated  by  the  integers  n = 1, 

2,  .... 


III.  NUMERICAL  PROCEDURE  FOR  SOLVING  EIGENVALUE  PROBLEM 

The  eigenvalue  problem  defined  by  the  sixth-order  system  (9),  (10), 
and  (11)  or  (12)  bears  some  similarity  to  the  Orr-Sommerfeld  equation. 

In  particular  the  coefficient  of  the  highest  derivative  contains  Re”1; 
in  liquid-filled  shell  applications  values  of  Re  ~ 106  are  commonplace. 
Our  numerical  procedure  uses  one  of  the  methods  developed  for  solving 
numerically  the  Orr-Sommerfeld  equation,  viz.,  orthonormalization. 

A.  Method  of  Superposition 

An  outline  of  the  method  will  be  given  after  a more  compact 
notation  is  introduced.  Matrices  will  be  denoted  by  upper  case,  under- 
scored symbols.  Express  (9) -(12)  in  matrix  form: 


Y = Y(r)  G 

03) 

Y (b)  E = 0 

(14) 

Y(l)  F = 0 

(15) 

where  Y is  the  solution  vector  Y = {y  ...,  y } a 1 x 6 matrix.  G is 

a 6 x 6 matrix;  E_  and  F_  are  6x3  matrices  the  elements  of  which  are 
obtained  from  (9) -(12);  b is  the  radial  coordinate  where  initial 
conditions  are  specified.  The  elements  of  E_  depend  on  the  choice  of 
boundary  conditions  (11)  or  (12) . Each  of  the  three  columns  of  E and  F 
has  only  one  non- zero  element. 
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The  general  solution  of  (13)  can  be  expressed  as  a linear  combina- 
tion of  six  linearly  independent  solutions.  If  these  are  chosen  so  that 
(14)  is  satisfied,  then  only  three  independent  solutions  are  available. 
These  are  denoted  by 

S.  (r)  = {s^,  ...  s^}  j = 1,  2,  3, 

so  that  the  solution  can  be  expressed  as 

Y(r)  = r_  S = YjSj  + y 2— 2 + Y3-3  (16) 


where  S is  a 3 x 6 matrix  whose  rows  are  the  vectors  , and  - 

(Y  Y , Y3^  is  a 1 x 3 matrix  of  constant  coefficients.  Note  that  (14) 

can  be  written  as 


Y (b)  E = _[_  S(b)  E = 0 


or 


S (b)  E = 0 


(17) 


The  eigenvalue  problem  is  solved  by  iteration.  A value  for  C is 
selected  at  each  stage  of  the  iteration  and  one  integration  pass  is  made 
between  r = b and  r = 1 for  each  of  the  three  independent  solutions, 

S_  , S_,  and  S_,  giving  the  s..(l),  j = 1,  2,  3 and  i = 1,  ....  6.  From 

(15) 


Y(l)  F = £ S(l)  F = 0 . (18) 

A non-trivial  solution  for  r_  requires 

D (C)  = Det  [S(l)  F ] = 0 . (19) 

The  iteration  process  provides  a systematic  method  of  choosing  the 
values  of  C used  in  the  numerical  integration  of  (9)  until  (19)  is 
satisfied  to  within  a specified  accuracy.  At  each  stage  of  the  itera- 
tion T is  obtained  by  solving  (18);  finally  the  eigenfunctions  are 
found  from  (16) . 

The  method  of  superposition  described  above  cannot  be  carried  out 
successfully  for  large  values  of  Re.  The  initially  linearly-independ- 
ent  solutions,  S^,  S^>  anc^  .8.3’  become  numerically  dependent  during  the 
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numerical  integration.  The  determinant  in  (19)  then  generally  fails  to 
yield  a clear-cut  root,  C,  and  the  iterative  process  does  not  converge. 
This  occurs  because  one  or  more  of  the  solutions  grows  rapidly  with  r, 
yielding  numbers  of  higher  orders  of  magnitude  than  the  other  solutions. 
A method  of  orthonormalization  is  employed  in  this  work  to  maintain 
linear  independence  and  reasonable  order-of-magnitude  solutions.  The 
technique  which  we  adopt  is  based  on  a version  of  Godunov’s  procedure7 
described  by  Conte8. 

B.  Numerical  Integration  and  Orthonormalization 


The  linear  independence  of  the  solutions  can  be  maintained  by  the 
orthogonalization  process  described  below.  The  solutions  can  be  further 
conditioned  by  periodically  renormalizing  them,  thus  preventing  any  of 
them  from  growing  inordinately  large.  Procedures  to  test  the  solution 
vectors  for  loss  of  independence  are  available9;  however,  to  avoid 
unduly  complicating  the  computer  program,  we  specify  the  number  and 
locations  of  orthonormalizations. 


The  integration  procedure  starts  with  orthonormal  initial  values 
(discussed  in  next  section)  specified  for  s..  at  r = 0.  When  b = 0, 

an  analytical  solution  near  r = 0 must  be  used  to  obtain  the  solution 
vectors  at  the  "matching  radius"  e;  when  b =£  0,  the  value  of  e is  taken 
equal  to  b. 


The  numerical  integration  of  (9)  is  carried  out  over  the  interval 
e < r < 1 using  a fourth-order  Runge-Kutta  technique.  This  interval  is 
divided  into  N subintervals  of  equal  length,  (1  - e)/N.  At  the  end  of 
each  subinterval  the  solution  vectors  are  orthonormal i zed  using  the 
Gram-Schmidt  process.  At  the  end  of  the  Xth  subinterval  there  are 


three  solution  vectors  S_^, 


and  S^,  where  X = 1, 


N denotes  the 


subinterval  index.  For  the  complex  solution  vectors  we  define  a non- 
standard inner  product 


6 


7.  S.  Godunov , "On  the  Numerical  Solution  of  Boundary -Value  Problems 
for  Systems  of  Linear  Ordinary  Differential  Equations ,"  Uspekhi 
Mat.  Nauk .,  Vol.  16,  1961,  pp.  171-174. 

8.  S . D.  Conte,  "The  Numerical  Solution  of  Linear  Boundary  Value 
Problems, " SIAM  Review,  Vol . 8,  1966,  pp.  309-321. 

9.  M.  R.  Scott  and  H.  A.  Watts,  "SUPORT-A  Computer  Code  for  Two-Point 
Boundary -Value  Problems  via  Orthonormalization, " Sandia  Labora- 
tories Report  SAND  75-0198,  Albuquerque,  New  Mexico,  June  1975, 
pp.  89-94. 
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The  two  vectors  are  orthogonal  when  • £2  = 0.  The  norm  of  a vector 
is  defined  by 


This  inner-product  was  used  by  Davey1®;  it  has  the  advantage  that  the 
solutions  are  analytic  functions  of  C,  a useful  property  in  the  itera- 
tion process.  It  has  the  disadvantage  that  the  norm  can  be  zero  for  a 
non-zero  vector;  however,  this  difficulty  has  never  arisen. 

The  new  set  of  orthogonal  vectors,  denoted  by  and  S^,  is 

obtained  using  the  following  transformation: 

Sx  = AX  Sx  (20) 


where  AX  is  a 3 x 3 matrix  whose  elements  are  given  by 

X _ X _ X _ . 

“ll  a22  a33  ' 1 


X _ X _ X _ _ 
a21  a31  a32  ° 


aX2  = QX  [(SX  • SX)(SX  • SX)  - (SX 


ij)  t§.j 


§.2>] 


(21) 


“i3  = QX  tfij  ‘ — 1^  1—2  ’ 


CiJ 


5.2)1 


a 


23 


CS.2  • & / csx  • Sx) 


where 


1/[(SX  • Sx)2 


<5i 


sb  (Sx 


■2J  '•-3  -3 


sx)] 


The  solution  vectors  are  then  normalized,  forming  three  orthonormal 
basis  vectors  T^  (j  = 1,  2,  3).  The  complete  orthonormalization  pro- 
cess  can  be  expressed  as 


10.  A.  Davey , "A  Simple  Numerical  Method  for  Solving  Orr-Sormerfeld 

Problems , " Quarterly  Journal  of  Mathematics  and  Applied  Mechanics , 
Vol.  26,  Part  4,  1973,  pp.  401-411. 
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(22) 


The  elements  of  TA  are  denoted  by  tA.; 

li 


the  elements  of  the  matrix  B are 


en  - 


si 


e22  = 


sA 

-2 


*33  = l/ 


si 


*21  = 4 = *32  = ° 


qA  _ X , 
*12  °12^ 


sA 

-1 


(23) 


„A  _ A . 
*13  a13^ 


sA 

-1 


flA  X / 
*23  a23/ 


sA 

-2 


The  three  solutions  T.  are  evaluated  at  the  end  of  the  Ath  sub- 


-1 


,A+1 


for  the  continuation  of 


interval  and  taken  as  initial  values  of 

the  numerical  integration  into  the  (A  + l)st  subinterval.  This  process 
of  integration  and  orthonormalization  is  carried  out  for  each  successive 

subinterval,  finally  yielding  three  orthonormal  bas_s  vectors  J^(l)  • 


C.  Initial  and  Terminal  Conditions 


The  initial  conditions  of  r = e must  be  selected  so  that  the 
boundary  conditions,  (11)  or  (12),  are  satisfied.  ?or  a cylinder  with 
an  inner  concentric  rod  the  initial  conditions  at  r = e are  taken  to  be 


SjCe)  = {0,0, 1,0, 0,0}  , S2(e)  = {0,0, 0,0, 1,0} 

S3(e)  = {0,0, 0,0, 0,1}  . 


(24) 


These  orthonormal  vectors  satisfy  (11),  or  equivalently,  (17). 

For  a cylinder  with  no  inner  rod  boundary  conditions  (12)  are 
specified  at  r = 0.  The  integration  of  (9)  must  be  carried  out  in  two 
stages  because  the  inverse  powers  of  r prevent  direct  numerical  inte- 
gration from  r = 0.  The  procedure  we  adopt  is  to  obtain  three  linearly 
independent  analytical  solutions  near  the  axis  which  satisfy  (12)  at 
r = 0.  Numerical  values  of  these  solutions  a~  the  matching  radius 
r = e are  used  as  initial  conditions  for  numerical  integration  over  the 
interval  e < r < 1.  The  analytical  solutions  are  obtained  from  power 
series  expansions  for  u,  v,  w,  and  p using  the  method  of  undetermined 
coefficients  to  obtain  recursion  formulas.  The  power  series  coeffi- 
cients for  m = 1 are  presented  in  Appendix  B . There  are  three  arbitrary 
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coefficients  (bQ,  d^,  e^) , which  are  specified  so  as  to  give  three 

independent  solutions.  Terminal  conditions  (10)  lead  to  (19),  modified 
by  orthonormalization  so  that  T(l)  replaces  £(1)  . In  terms  of  T(l)  the 
following  equation  must  be  solved  for  C: 


*11 0) 

tNnCl) 

‘si'1) 

D(C)  = 

t^Cl) 

*&(« 

= o . 

(25) 

‘>1 

*24^ 

*>> 

Once  a solution  to  (25)  has  been  found,  y^,  y^,  and  y^  of  (16)  can 

be  determined  to  within  a constant  multiple.  The  matrix  from  which  the 
determinant  in  (25)  originates  is  generally  of  rank  2.  Thus,  we  set 
y^=l  so  that  Yj  and  y^  are  found  from 

tll ('1)  Y1  + t21  1 Y2  ^31  1 

(26) 

tl2(1)  Yj  + t22^1^*  y2  t32('1^  * 


D.  Iteration  Procedure  to  Solve  for  Eigenvalue 

An  eigenvalue  C = + iC^  is  a root  of  (25) , which  is  equivalent 

to  the  following  system  of  two  real  equations: 


WP  = 0 > Wi3  = 0 


(27) 


where  and  are  real  and  imaginary  parts  of  D.  This  system  is 
solved  by  a method  of  successive  iteration.  An  initial  guess  is 

assumed  which  must  be  reasonably  close  to  the  solution  because  of  the 

sensitivity  of  D and  D to  variations  in  Cn  and  CT . An  iteration 
K 1 K 1 

operator,  J,  is  applied  to  yielding  a new  estimate  C^,  etc.  A 
sequence  of  iterates  is  determined  from 

C , = J (C  ) . 

p+1  v \iJ 
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When  (25)  is  satisfied  and  |c^+j.  - C | is  smalL  to  within  specified 

tolerances,  for  some  y,  the  iteration  is  stopped.  We  found  it  advisable 
to  satisfy  both  requirements  since  the  condition  on  C could  sometimes  be 
satisfied  without  the  one  on  D being  satisfied. 


The  iteration  operation  used  here  is  an  extension  of  Newton's 
method.  For  small  changes  in  C and  D between  the  yth  and  (y+l)st 
iterations  we  can  write 


R - Dr  = (9V9Cr}  [cr  - CR  ] + CSD^/aCj)  [Cj  - Cj  ] 


y+1  y 


y+i  y 


y+1  y 


(28) 


DI  , " DI  = C3Dj/3Cr)  [Cr  - CR  ] + ODt/3Ct)  [Ct  - CT  ] 


y+1  y 


y+1  y 


I'  r 1 i , i 
y+l  y 


The  partial  derivatives  are  approximated  by 


3V3CR  - (DR  <CR  * iCR’CI  > - °R  <CR  -CI 
y y y y 


(29) 


aDj^/ac-,.  = [dr  (cr  ,Cj  + ac:)  - dr  (cr  ,Cj  )]/acj  , 
y y y y 


with  similar  expressions  for  8DT/oCr>  and  9DT/9C  . Appropriate  small 

IK  11 

values  for  ACD  and  ACT  are  used,  e.g.  ACD  = SL  C unless  Cn  is  small,  in 
K 1 K K K 

which  case  CR  = £;  i is  an  input  parameter,  typically  SL.=  0.001.  The 

goal  of  each  iteration  is  to  make  D , vanish.  We  set  D„  = DT  = 0 

y+1  R . I , 

y+1  y+1 

in  (28)  and  thereby  get  formulas  for  C and  CT  in  terms  of  C 

y+1  y+1  y 

and  Cj  . 
y 

A simplification  is  possible  at  this  point,  since,  with  the  chosen 
definition  of  inner  product,  D(C)  is  an  analytic  function  of  C which 
satisfies  the  Cauchy- Riemann  conditions: 


9Dr/9Cr  = 9DJ/9CJ  , 9V9CI  = " 9DI/9CR  * ^ 

We  substitute  (30)  into  (28)  (with  D)J+^  = 0)  &nd  obtain  the  operator 
J(C^)  as  the  solution  to  the  set  of  equations: 
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(31) 


(9VaV  CR  , - '3  W CI  , - H!  fVV 

y+1  y+1 

(3D  /3C  ) CR  ♦ (3Dr/3Cr)  C - H2  (CD) 

y+1  y+1 

where 

H1  <VV  s (3DR/3CR)  Cr  - (3Di/3Cr)  C, 

y y 

H2  (CD)  = (3D  /3C  ) CR  * (3Dr/3Cr)  C, 

y y 

Since  (31)  and  (32)  contain  only  partial  derivatives  with  respect 
to  CR,  DR  (Cr  , Cj  + ACj)  in  (29)  and  (CR  , Cj  + ACj)  in  the 

y y y y 

expression  for  9Dj/9Cj  need  not  be  evaluated.  Practically,  this  cuts 

the  required  computer  time  almost  in  half  since  only  six  integration 
passes  (two  for  each  of  the  three  independent  solutions)  are  needed  per 
iteration,  versus  the  twelve  that  would  be  required  if  D(C)  were  not 
analytic . 

E . Evaluation  of  Eigenfunctions 

For  some  applications  we  may  also  wish  to  evaluate  the  flow 

variable  perturbations  u(r),  v(r),  w(r),  and  p(r),  or,  equivalently, 

Y(r)  = {y. , . ..,  yA>  the  eigenfunctions  associated  with  the  eigen- 
— Id 

value  C.  The  solution  Y_  is  expressed  as  a combination  of  the  three 
linearly  independent  solutions,  according  to  (16).  Because  of  the 
orthonormalization  at  the  end  of  each  subinterval  the  constants  in 

(16)  are  different  for  each  sub interval:  yj,  y^,  Yj*  The  solutions 

S^(r),  and  S^CO  are  comPuted  over  each  subinterval  and  stored; 

the  evaluation  of  the  eigenfunctions  reduces  to  the  problem  of  deter- 

A A , A 
mining  y , y , and  y . 

At  the  (A  - l)st  junction  point  r = e + (1  - e) (A  - 1)/N  and 
the  solution  is 

YX_1  (rx_p  - rX_1  SX_1  (rx  i)  . (33) 


- D 


R 


- 


(32) 
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The  initial  value  of  the  solution  in  the  Xth  subinterval  is  (see  Section 
III .B) 


-X  <rx-d  = - l*"1  ^X-l3 


At  each  junction  point  = T' ^ ^ and 


sx  = ba  1 sx_1 


(34) 


from  (22).  The  elements  of  the  matrix  B^  \ given  by  (23),  are  stored 
for  all  subintervals  during  each  iteration.  Continuity  requires  that 

YA_1  (r.  ,)  = YA  (r^  ).  Equating  the  right-hand  sides  of  (33)  and  (34) 

gives 


r*-1  s»-i 


* - 


Bx-x  sx_1 


or 


rA  1 = rx  ba_1 


(35) 


which  is  used  to  calculate  the  T A ^ from  known  T_  . At  the  last  (Nth) 

subinterval  y!?,  and  are  known.  In  our  procedure  y^  = 1,  and 

y^  and  y ^ are  obtained  from  (26).  The  y^  (j  = 1,2,3)  are  found 

using  (35),  and  so  forth,  down  to  X « 2.  The  appropriate  linear  combin- 
ations of  S^,  S^,  and  are  evaluated  in  each  sub interval  to  give  the 

eigenfunctions  over  the  interval  e < r < 1 . 


F . Additional  Aspects  of  the  Computation 

As  mentioned  in  Section  III.D,  convergence  of  the  iteration  pro- 
cess requires  that  the  initial  guess  for  the  eigenvalue  C be  reasonably 
close  to  the  actual  value,  say  within  10%;  typically,  convergence  cri- 
teria are  satisfied  in  five  or  less  iterations.  For  solid  body  rotation 
and  m = 1,  satisfactory  trial  values  can  be  obtained  from  Wedemeyer’s 
simple  formulae  for  eigenvalues1 * 3 obtained  by  applying  a viscous 
correction  (described  in  Section  IV)  to  Stewartson’s  results2.  In  fact, 
the  latter  can  also  be  used  for  the  initial  gusss,  with  = 0 because 

the  perturbation  is  inviscid.  The  are  given  in  tabular  form  in 

Reference  1 for  given  c and  k and  radial  mode  numbers  n = 1,  2,  3.  (In 
Reference  1,  CD  is  denoted  by  t , and  2j  = k - 1.)  The  value  of  n is 

the  maximum  number  of  (non-boundary)  zeros  in  the  real  and  imaginary 
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parts  of  the  set  of  u,  v,  w,  and  p.  A typical  solution  is  illustrated 
in  Figure  1*  which  shows  the  real  part  of  p(r)  for  the  lowest  two 
radial  modes  for  the  case  c = 3.148,  Re  = 8977,  k = 3,  m = 1.  The 
number  of  zeros  is  n - 1 here,  but  for  v the  number  of  zeros  is  n. 

In  Figure  2 the  real  part  of  w is  plotted  to  demonstrate  the 
presence  of  the  boundary  layer  on  the  sidewall. 

It  was  stated  in  Section  III.C  that  three  independent  analytic 
solutions  to  (9)  were  obtained  near  r = 0.  Though  it  has  not  been 
proven,  it  is  fairly  certain  that  any  other  independent  solution  would 
be  singular  at  r = 0,  and  thus  inadmissable.  For  m = 1 we  obtain  three 
independent  solutions  by  choosing  the  complex  constants  bQ,  d^,  and  e^ 

in  Appendix  B in  the  following  way:  during  each  iteration  {bQ,  d^,  e^} 

is  equal  to  {1,0,0}  for  the  first  integration  pass,  {0,1,0}  for  the 
second  pass,  and  {0,0,1}  for  the  third.  Approximately  ten  terms  of  the 
series  are  used  to  approximate  the  functions  for  e = 0.001.  The  actual 
number  of  terms  used  is  adjusted  to  insure  that  the  series  converge 
with  a difference  between  successive  partial  sums  of  no  greater  than 
typically  1 x 10-6. 

Since  m = 1 modes  are  the  only  ones  which  lead  to  an  overturning 
moment  on  the  projectile,  only  these  will  be  considered  in  this  report. 
However,  other  modes,  particularly  m = 0,  have  also  been  computed. 


IV.  ENDWALL  BOUNDARY  CONDITION  CORRECTION 

The  technique  described  in  Sections  II  and  III  will  provide 
natural  frequencies  within  the  limitations  of  the  linearized  theory. 
However,  as  pointed  out  in  Section  II,  in  order  to  perform  a modal 
analysis  in  the  longitudinal  direction  it  is  necessary  to  abandon  the 
no-slip  boundary  conditions  u = v = 0 on  the  endwalls;  only  w = 0 
could  be  satisfied.  The  perturbation  velocities  thus  satisfy  the 
inviscid  boundary  conditions.  Because  these  incorrect  boundary 
conditions  have  a significant  effect  on  C^  a correction  was  devised, 

analogous  to  Wedemeyer's  viscous  correction1*3  to  Stewartson's  inviscid 
results.  This  correction  is  briefly  described  here;  details  are  given 
in  Appendix  A. 

Experiments  with  a liquid-filled  gyrostat  revealed  discrepancies 
between  observations  and  Stewartson's  inviscid  analysis.  To  account 
for  these  Wedemeyer  analyzed  the  perturbation  boundary  layers  along 
the  sidewall  and  endwalls.  He  concluded  that  the  viscous  flow  problem 
is  equivalent  to  an  inviscid  problem  with  boundary  conditions 


* Figures  1 and  2 are  located  on  page  26. 
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at 


(36) 


* 

u 


0 


r = 1 - 6a 

w 


# = 0 at 


z = 6c  and 

w 


£ = (2c  - 6c  ) , 

w 


where  6a  and  6c  are  non-dimensional  displacement  thicknesses  for  the 
w w 

sidewall  and  endwall  boundary  layers,  respectively.  The  complex  quanti- 
ties 6a  and  6c  are  algebraic  functions  of  Re  and  the  inviscid  CL  but 
w w & R 

are  independent  of  r,  0,  and  t.  The  analysis  is  applicable  for  large 
Re  (but  not  so  large  as  to  cause  boundary  layer  instabilities) . Large 
Re  requires  that  l^a^]  <<  1 and  1 6 | <<  1,  conditions  which  are  gener- 
ally true  for  Re  > 103.  Solving  the  inviscid  problem  with  (36)  modi- 
fies the  eigenfrequency  by  a small  complex  increment.  The  resulting 
frequency  and  damping  rate  agree  with  experimental  observations. 

In  our  analysis  the  viscous  diffusion  terms  are  present  in  the 
perturbation  equations.  The  correction  is  needed  here  because  condi- 
tions u = v = 0 on  z = 0,  2c  are  not  satisfied.  This  correction,  in 
principle,  is  different  from  Wedemeyer’s;  but  after  proper  formulation 
the  analysis  proceeds  in  much  the  same  way.  A sidewall  correction  is 
not  required  since  the  perturbation  solution  satisfies  the  no-slip 
conditions  at  r = 1.  A boundary  layer  analyses  is  made  on  the  differ- 
ence between  the  flows  with  and  without  the  no-slip  conditions  on  the 
endwalls,  and  a displacement  thickness  6c^  is  obtained  which  has  the 

same  properties  as  6c^.  When  = 0,  6c^  = 6c^.  The  eigenvalue 

problem  is  solved  again,  now  using  the  modified  complex  half-height, 

c - 6c  . No  computational  difficulties  are  introduced  because  6c 
s s 

directly  affects  only  K in  (7) , which  appears  in  the  already  complex 

(9) . The  effects  of  this  correction  on  C and  C are  discussed  in 

Section  V. 


V.  DESCRIPTION  OF  NUMERICAL  RESULTS 

A program  as  extensive  as  the  eigenvalue  computation  program 
requires  checks  to  establish  its  reliability.  Three  types  of  checks 
will  be  treated:  (a)  related  problems  for  which  results  are  available 
in  the  literature;  (b)  Wedemeyer!s  approximate  solution  to  the  eigen- 
value problem;  (c)  experimental  data. 

A.  Check  With  Flow  Between  Rotating  Cylinders 


Many  investigations  have  been  carried  out  on  the  stability  of 
flow  between  concentric  rotating  cylinders,  often  called  Taylor’s 
problem;  when  the  flow  is  unstable  Taylor  vortices  are  generated.  The 
basic  flow  is  steady  but  not  solid-body  rotation  because  the  angular 


20 


velocities  of  the  two  cylinders  differ.  V(r)  has  r and  r“  terms,  and 
is  independent  of  z for  infinite  length  cylinders.  Since  our  program 
is  written  for  general  V(r),  eigenvalues  can  be  calculated  for  the  two- 
cylinder  problem. 

The  first  comparison  is  made  with  a result  in  Reference  11  for 
which  the  parameters  are  Re  = 796.8,  b = 0.70,  £2^/^a  = -1,  K = 13.280, 

and  m = 2,  where  b is  the  non-dimensional  radius  of  the  inner  cylinder, 
and  ft,  and  ft  are  the  dimensional  angular  velocities  of  the  inner  and 

outer  cylinders,  respectively.  The  computed  eigenvalue  is  CR  = 0.6763 

and  Cj  = -0.0007,  compared  to  CR  = 0.6759  and  Cj  = 0.0000  in  Reference 

11.  The  second  comparison  is  made  with  a result  of  Reference  12  for 

Re  = 4270.2,  b = .9512,  «,/fi  = 0,  K = 64.10,  and  m = -1.  Values  of 

d a 

C = -0.5141  and  C = -0.1443  are  obtained  both  from  our  program  and 
R I 

from  Reference  12.  The  agreement  between  our  results  and  those  of 
References  11  and  12  lends  support  to  the  validity  of  all  of  them. 

B.  Check  With  Wedemeyer's  Solution  and  Parametric  Study  of  Eigenvalues 

No  previous  viscous  perturbation  calculations  exist,  but  Wedemeyer’s 
viscous  correction  to  the  Stewartson  theory  provides  an  approximate 
solution  to  the  eigenvalue  problem  which  agrees  well  with  experiment. 
Therefore,  results  from  his  theory  will  be  compared  with  ours;  a suffi- 
cient number  of  cases  will  be  presented  to  exhibit  the  variation  of  the 
eigenvalue  with  Re  and  c. 

Tables  1-3  provide  the  eigenvalues  for  the  m ■ 1,  k = 3,  n = 1 
mode  over  the  ranges  c = 3,  4,  5 and  Re  = 10^,  10^,  10^,  and  10^.  The 
endwall  displacement  thicknesses  6c^  and  are  also  presented;  there 

is  no  reason  for  these  two  quantities  to  be  the  same;  but  in  fact  they 
differ  at  most  by  a few  percent.  The  effect  of  the  boundary  condition 
correction  on  CD  and  CT  is  demonstrated  in  these  tables.  Except  for 

c = 3 and  the  two  smaller  Re,  the  percentage  correction  to  Cj  is  much 

larger  than  that  to  CR.  With  boundary  condition  correction  the  CR  differ 

from  Wedemeyer's  results  by  a few  percent  at  most  for  all  cases.  How- 
ever, for  the  three  aspect  ratios  and  Re  = 10^  and  101*  the  differences 

11.  E.  R.  Krueger,  A.  Gross , and  R.  C.  DiPrima , "On  the  Relative 
Importance  of  Taylor-Vortex  and  Non-Axisymmetric  Modes  in  Flow 
Between  Rotating  Cylinders , " J_.  Fluid  Mech. , (1966) , Vol.  24, 

Part  3,  pp.  521-538. 

12.  P.  M.  Eagles,  "On  Stability  of  Taylor  Vortices  by  Fifth-Order 
Amplitude  Expansions , " J . Fluid  Mech. , (1971) , Vol.  49,  Part  3, 
pp.  529-550. 


in  Cj  are  greater  than  10%,  and  for  c = 5,  Re  = 103  it  is  43%.  In 

principle,  our  results  should  be  more  accurate  since  no  correction  is 
needed  at  the  sidewall  and  viscous  diffusion  effects  in  the  interior  are 
included.  Whether  or  not  this  is  the  case  requires  a comparison  with 
experiment  (or  perhaps  a full  numerical  solution  to  the  perturbation 
problem) . On  the  basis  of  the  following  discussion  it  is  concluded  that 
sufficient  experimental  evidence  to  discriminate  between  Wedemeyer's 
results  and  ours  does  not  exist. 


Table  1.  Eigenvalues  and  Half-Height  Corrections  for  Solid-Body  Rotation: 
m = 1,  k = 3,  n = 1,  c « 3, 


Re 

r 

B.C. 

Correction 

No  B.C. 
Correction 

f 

B.C. 

Correction 

No  B.C. 
Correction 

= a2ft/v 

6c 

w 

6c 

s 

^ed. 

CR 

CR 

UI 

Wed. 

ci 

CI 

106 

.000855- 

.000856- 

.00472 

.00484 

.00509 

.000996 

.00100 

.000630 

.00127  i 

.00127  i 

10s 

.00270- 

.00271- 

.00553 

.00564 

.00643 

.00315 

.00324 

.00207 

.00401  i 

.00401  i 

104 

.00855- 

.00858- 

.00808 

.00815 

.01066 

.00996 

.01100 

.00727 

.01267  i 

.01273  i 

103 

.02705- 

.02754- 

.01614 

.01578 

.02382 

.03150 

.04256 

.03052  ’ 

.04008  i 

.04069  i 

Table  2.  Eigenvalues  and  Half-Height  Corrections  for  Solid-Body  Rotation: 
m * 1,  k * 3,  n * 1,  c 3 4. 


B.C. 

No  B.C. 

B.C. 

No  B.C. 

Re 

£ 

Correction 

Correction 

c 

Correction 

Correction 

= a2fl/v 

6c 

w 

6c 

s 1 

'Sved. 

CR 

CR 

Cj 

wed. 

ci 

ci 

.000805- 

.000806- 

i 

106 

.00150  i 

.00151  i 

.24352 

.24375 

.24390 

.000915 

.000927 

.000636 

.00255- 

.00256- 

10s 

.004  76  i 

.00477  i 

.24454 

.24475 

.24525 

.00289 

.00300 

.00208 

10“ 

.00805- 

.00815- 

.24775 

.24797 

.24953 

.00915 

.01024 

.00731 

.01505  i 

.01515  i 

.02547- 

.02685- 

103 

.04758  i 

.04859  i j 

.25793 

.25772 

.26277 

.02894 

.03995 

.03051 
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Table  3.  Eigenvalues  and  Half-Height  Corrections  for  Solid-Rody  Rotation: 
m = 1,  k=  3,  n = 1,  c = 5. 


B.C. 

No  B.C. 

B.C. 

No  B.C. 

Re 

Correction 

Correction 

ci 

‘Wed. 

Correction 

Correction 

= a2n/v 

6c 

w 

6c 

s 

'Sved. 

CR 

CR 

ci 

ci 

.000783- 

.000784- 

106 

.00181  i 

.00181  i 

.40225 

.40227 

.40237 

.000822 

.000837 

.000602 

.00248- 

.00249- 

105 

.00573  i 

.00574  i 

.40330 

.40333 

.40364 

.00260 

.00272 

.00198 

.00783- 

.00800- 

104 

.01811  i 

.01827  i 

.40664 

.40666 

.40768 

.00822 

.00935 

.00699 

.02476- 

.02723- 

103 

.05727  i 

.05875  i 

.41719 

.41674 

.42010 

.02600 

.03729 

.02966 

C.  Check  With  Experimental  Data 

This  check  is  the  most  useful  one.  Unfortunately  there  are  not 
sufficient  data  for  completely  filled  cylinders  to  make  the  necessary 
crucial  comparisons.  Most  of  the  experimental  data  were  generated  in 
a gyrostat;  the  technique  of  determining  C and  CT  is  discussed  in 

Reference  1.  The  were  not  measured  directly  but  were  determined 

from  the  observed  maximum  rate  of  divergence  of  the  gyrostat  motion  and 
relationships  taken  from  Wedemeyer's  theory;  the  latter  fact  indicates 
a possible  logical  inconsistency  in  a comparison  of  Wedemeyer's  theory 
with  the  experimental  C^ . 

Results  for  two  values  of  Re  are  compared  in  Table  4 for  k = 3, 
m = 1,  n = 1,  and  c = 3.148;  the  data  in  Case  I come  from  Reference  13; 
Case  II  was  done  specifically  for  this  report.  Both  Wedemeyer's  and 
our  results  for  CD  agree  with  the  experimental  values  for  both  cases  to 

within  the  estimated  uncertainty.  The  data  were  not  accurate  enough  to 
obtain  C^  for  Case  I,  Re  = 5.2  x 105;  but  from  results  for  partially- 

filled  cylinders  we  should  expect  agreement  at  this  Re.  For  Case  II 
the  theoretical  and  experimental  results  for  C^  agree,  considering  the 

accuracy  of  the  latter.  Thus,  it  is  an  open  question  whether  our 
results  are  in  fact  superior  to  Wedemeyer's. 


13.  W.  P.  D^Amioo^  Jr.  j private  oorrmunieation. 
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Table  4.  Comparison  of  Predicted  and  Measured  Eigenfrequencies 
for  Mode  m = 1,  k = 3,  n = 1 for  Basic  Flow  in  Solid- 
Body  Rotation. 


Case  I : Re 

= 5.195  x 105 

Case  II: 

Re 

= 8.977  x 103 

CR 

CI 

CR 

CI 

Gyrostat  Measurements 

0.048  ± 

(Insufficient 

0.051  ± 

0.011  ± 

.002 

Accuracy) 

.003 

.002 

Stewartson  Inviscid 
Perturbation 

0.0462 

0 

0.0462 

0 

Wedemeyer  Correction 
to  Stewartson  Value 

0.0468 

0.0014 

0.0500 

0.0104 

Result  from  Present 
Viscous  Perturba- 
tion— No  Correction 

0.0473 

0.0009  j 

0.0530 

0.0078 

Result  from  Present 
Viscous  Perturba- 
tion --Correct  ion  5c 

s 

0.0469 

0.0014  , 

0.0505 

0.0116 

Most  of  the  experiments  were  conducted  with  partially-filled 
cylinders  because  these  are  the  more  practical  cases  and  the  filled 
cylinder  experiment  is  more  difficult.  Wedemeyer3  compared  his  calcu- 
lations with  data  obtained  from  the  gyrostat  experiments,  and  the  agree- 
ment was  impressive.  The  data  were  obtained  for  a cylinder  of  nominally 
85%  fill-ratio,  but  the  actual  value  was  varied  in  the  experiments.  The 
viscous  correction  theory  of  Wedemeyer  requires  100%  fill -ratio,  as  does 
ours.  But  since  he  applies  his  correction  to  Stewartson *s  results  for 
partially-filled  cylinders,  he  evidently  assumes  that  it  is  applicable 
there  also;  this  assumption  seems  justified  by  the  favorable  comparison. 
There  is  no  counterpart  to  Stewartson* s results  to  which  our  boundary 
condition  correction  can  be  applied;  thus  a clear-cut  check  of  the  two 
theoretical  results  against  experimental  data  is  not  possible. 

Analysis  of  the  available  data  and  calculations  indicates  that  the 
difference  between  the  C^*s  for  85%  and  100%  ratios  is  negligible  for 

Re  ^ 10^  and  c = 3;  however,  for  Re  ^ 101*  the  effect  of  fill-ratio  on 
Cj  is  more  pronounced.  For  example,  with  Re  = 5.05  x 102,  c=  3.08, 

k = 3,  m = 1,  and  n = 1,  Wedemeyer* s formulae  give  = 0.0490  (85%) 

and  0.0437  (100%).  Whether  or  not  these  accurately  represent  the  effect 
of  fill-ratio  at  low  Re  is  not  known,  and  appropriate  experimental  data 
are  not  available.  Further  experimentation  is  required. 
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VI.  CONCLUSIONS 


Although  the  main  objective  of  the  work  reported  here  is  to  develop 
methods  to  compute  the  eigenvalues  for  small  perturbations  of  spin-up 
flow  in  a cylinder,  this  Part  I is  restricted  to  consideration  of  solid- 
body  rotation  so  that  the  numerical  technique  can  be  described  in  detail 
without  the  distraction  of  the  complexities  that  arise  in  spin-up,  and 
yet  provide  results  that  can  be  checked  with  other  theoretical  work  and 
experimental  data.  The  perturbation  equations  are  solved  using  an 
orthonormalization  technique  to  insure  linear  independence  of  the 
solutions,  and  the  eigenvalues  are  determined  iteratively.  A correction 
for  the  endwall  boundary  condition  is  derived  and  incorporated  into  the 
program  which  then  provides  eigenvalues  that  agree  with  other  results 
in  the  literature  and  with  experimental  data. 
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n = 1 


Figure  1.  Real  part  of  p(r)  for  two  radial  modes,  with  c - 3.15, 
Re  = 8977,  m = 1 , k = 3. 


Figure  2.  Real  part  of  w(r)  for  n = ',  with  c = 3.15, 
Re  = 8977,  m = 1 , k = 3. 
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LIST  OF  SYMBOLS 


Length,  velocity,  pressure,  and  time  are  non-dimensionalized  by  a, 
a ft,  pft2a  , and  ft-1,  respectively. 

a cross-sectional  radius  of  cylinder  [m] 

3x3  orthogonalization  matrix  in  Ath  sub interval  (see  (20)) 
b value  of  r where  initial  conditions  are  specified 

3x3  orthonormalization  matrix  in  Ath  subinterval  (see  (22)) 

c half-height  of  cylinder 

C E CD  + iCT,  complex  eigenvalue 

K 1 

Cj  disturbance  decay  rate/ft 

CD  disturbance  frequency/ft 

K 

D(C)  e D_  + iDT,  3x3  determinant  whose  zeros  are  the  eigenvalues 
R 1 

(see  (19)) 

E,F  6x3  matrices  specifying  initial  and  terminal  conditions, 

respectively  (see  (14)  and  (15)) 

F symbol  denoting  a solution  u,v,w,p  (Appendix  A) 

(I  6x6  matrix  of  coefficients  in  (9)  (see  (13)) 

J iteration  operator  --  C ^ = J(C^)  (see  (31)) 

k,m,n  axial,  azimuthal,  and  radial  indices,  respectively,  of  a 
disturbance  mode 

K = kir/  (2c)  (see  (7)) 

N number  of  subintervals  in  orthonormalization 

p pressure 

pressure  perturbation  (see  (1)) 
radial  part  of  pressure  perturbation  (see  (6)) 


29 


LIST  OF  SYMBOLS  (Continued) 


p 

unperturbed  pressure 

r 

radial  coordinate 

Re 

= fta2/v  , Reynolds  number 

Sjl*  Sj6 

elements  o£  solution  vector  S. 

—j 

S. 

~J 

eigenvalue  problem  solution  vector,  j = 1,  2,  3 
(see  (16)) 

s 

3x6  solution  matrix  whose  rows  are  S_.  (see  (16)) 

orthogonalized  solution  matrix  in  Ath  subinterval 
(see  (20)) 

t 

time  [non-dimensional ] 

A 

t.  . 

31 

elements  of  3x6  matrix 

_A 

T 

orthonormal i zed  solution  matrix  in  Ath  subinterval 
(see  (22)) 

u,  v,w 

radial,  azimuthal,  and  axial  components,  respectively, 
of  velocity 

*/*  * 

/ 

ll,V,W 

radial,  azimuthal,  and  axial  components,  respectively, 
of  perturbation  velocity  (see  (1)) 

U,  V,  w 

radial  parts  of  u,  v,  and  w,  respectively  (see  (6)) 

U,  V,  W 

radial,  azimuthal,  and  axial  components,  respectively, 
of  unperturbed  velocity 

y 

boundary  layer  normal  coordinate  at  endwall  (see 

A. 3a-d) 

^1 3 * * * J ^6 

dependent  variables  in  canonical  form  of  eigenvalue 
differential  equation  system  (see  (8)) 

Y 

solution  vector  (y^ , y^} 
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LIST  OF  SYMBOLS  (Continued) 


z 


X 

a.  . 
il 

V V Y3 

6c  , 6c 
s w 


A 

r 


e 

0 

V 

p 

a 


axial  coordinate  (=  0 at  base  of  cylinder) 
elements  of  (see  (21)) 

elements  of  £ (see  (23)) 
constant  coefficients  in  (16) 

cylinder  half-height  correction  due  to  endwall  boundary 
layer--present  theory,  Wedemeyer  theory,  respectively 

displacement  thickness  of  endwall  boundary  layer  (see 
A. 6) 

1x3  matrix  {y^,  Y2»  Y^) 

value  of  r where  numerical  integration  is  initiated 

azimuthal  angle  [radians] 

kinematic  viscosity  of  fluid  [m2/s] 

density  of  fluid  [kg/m3] 

spin  rate  of  cylinder  [radians/s] 


Superscript 

X sub interval  index 

~ endwall  boundary  layer  approximation 

' d/dr  in  (8)  and  (9) 


Subscript 

I 

R 

X 


imaginary  part  of  complex  number 
real  part  of  complex  number 
sub interval  index 
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LIST  OF  SYMBOLS  (Continued) 


y iteration  index  in  solution  of  D(C)  =0  (see  (28)) 

0 solution  which  satisfies  inviscic  boundary  conditions 
at  endwalls  (Appendix  A) 

1 correction  to  o solution  (Appendix  A) 

00  endwall  value  of  o solution 


Miscellaneous 

• inner  product  --  S_^  • S^ 

||  ||  norm  of  a vector  --  | j S 


s . s 
li  2i 


6 
Z 

i=l 

2 

= S • S 
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APPENDIX  A:  CORRECTION  FOR  THE  ENDWALL  BOUNDARY  CONDITIONS 


The  derivation  of  the  correction  for  endwall  boundary  conditions, 
discussed  in  Section  IV,  is  presented  here;  this  correction  is  needed 
to  satisfy  u=v=0onz=0,2c.  The  perturbation  equations  for  solid 
body  rotation  are  (3),  the  linearized  Navier-Stokes  equations: 

4*  4*  4* 

(A. la) 
(A.  lb) 
(A.lc) 
(A.  Id) 

L U L 

where 


r 

e 

z 

* 

* 

* 

* 

ut 

i 

CD 

+ 

2v 

= - p 

rr 

* 

* 

* 

V- 

+ + 

2u 

= - r"1 

t 

e 

* 

* 

* 

wt  + 

w0 

= - p 
rz 

- 92/3r2  + 

r” 

1 3/3r  + 

1 o * _ 9 * -0  * 

} “ -1-  CK 7^-  n t*  “ ^ 


± * 


.-2 


1 u2 


^ ^ 

The  boundary  conditions  are  u = v = w = 0atz  =*0,*2c  and  at  r = a. 
The  solution  to  the  problem  is  denoted  by  F = (u,  v,  w,  p) . 

As  discussed  in  Section  II  we  solve  the  same  equations 


u +u  -2v=-p  + Re-1  (V2  u - r-2  u - 2r-2  v ) (A. 2a) 

0 0 o r o v 0 o o' 

t 0 r 0 


v + v + 2u  = - r-1  p + Re"1  (V2  v - r"zv  + 2 r“z  u ) (A 

°t  °e  0 °e  o °e 


w + w = - p + Re-1  V2  w 

°t  °e 


fru)  + v +rw  =0, 
or  o„  o 

0 z 


(A. 2c) 

(A. 2d) 


but  with  boundary  conditions  u =v  =w  = 0 at  r = a,  and  w = 0 at 
J ooo  o 

z = 0,  2c.  but  u # 0 and  v # 0 there.  Let  F denote  the  solution  to 
J ’ o o o 

this  problem,  for  which  a modal  analysis  in  the  z direction  is  possible. 

Since  erroneous  values  for  Cj  were  obtained  with  these  incorrect 


.2b) 
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(inviscid)  boundary  conditions  and  a modal  analysis  for  F could  not  be 
made,  a correction  to  was  sought  which  satisfies  the  no-slip  boundary 

conditions  at  z = 0,2c. 

Subscript  1 is  used  to  denote  the  correction.  Since  F and  Fq 

differ  only  in  their  boundary  conditions  at  z = 0,2c,  we  expect  F - F^ 

to  differ  from  zero  only  near  the  end-walls.  The  correction  F^  should 

have  the  property  that  u+u  = v,  + v =w+w  = 0 at  z = 0,2c,  and 

lololo 

F^  0 or  F^  + F^  •+  F,  away  from  this  boundary.  Thus,  F^  = F - Fq  must 

have  a boundary  layer  character  near  z = 0,2c;  and,  because  of  the 
similarity  of  this  problem  to  Wedemeyer's,  the  same  Re”  2 scaling  is 
appropriate  to  this  boundary  layer.  The  boundary  layer  approximation  to 

F^  is  denoted  by  F^  and  is  an  approximation  to  the  desired  correction, 

but  only  in  the  boundary  layer;  elsewhere  the  correction  is  negligible. 

Under  the  boundary  layer  approximation  F^  = F - Fq  becomes  F^  = 

F - F^  where  F^  is  the  boundary  layer  approximation  to  the  solution  of 

(A.2a-d);  the  F,  obtained  from  this  relation,  is  the  corrected  F or  the 

o 

desired  F,  to  within  the  boundary  layer  approximation.  Let  6 be  a 

1 _i^ 

measure  of  the  thickness  of  the  boundary  layer  so  that  6^  = 0 (Re  2)  . 

Using  the  form  of  the  solution  for  Fq  in  (6a-d)  it  can  be  easily  shown 

that  Fq  satisfies  the  inviscid  equations;  i.e.,  (A.2a-d)  with  v = 0, 

with  an  error  of  OfRe”1)  except  near  r = a,  where  the  sidewall  boundary 
layer  induces  a viscous  contribution.  Also,  the  changes  in  u and  v 
across  6^  are 


Auq  = 0 (Re”1)  , Avq  = 0 (Re  1 ) 


Thus, 


u 


u 


ui = 


u 


o°° 


+ U- 


V 


V + V, 

o 1 


V + 
o00 


v., 


with  the  same  error  0 (Re-1),  where  sub  00  denotes  evaluation  on  the  end- 
walls  z = 0,2c.  The  conclusions  of  this  paragraph  are  needed  when  the 
displacement  surface  is  determined. 

The  boundary  layer  approximation  applied  to  the  equations  for 
F^  = F - F , yields 
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(A. 3a) 


ui,  * v - 2vi  * pi  ■ Re'1  ui 

t 0 r yy 


v1  + v + 2u1  + r"1  p = Re-1  v 
lt  i0  1 i9  iyy 


(A. 3b) 


Pi  = ° 

y 


(A. 3c) 


(ru  ) + v.  - (rw.)  = 0 

1 r 1.  ly 


(A. 3d) 


for  the  endwall  at  z - 2c,  where  y = 2c  - z is  the  boundary  layer  normal 
coordinate.  The  boundary  conditions  are 


u = - u 
1 o°° 


vl  = - vn„  > at  y = 0 , u = V = 0 at  y = 


w^  = 0 


Since  the  flow  is  generated  by  changing  from  the  no-slip  to  the 

inviscid  type  boundary  conditions  there  is  no  mechanism  to  cause 
pressure  gradients.  Thus,  p][  = 0 and  p = 0.  Then  (A.3a-b)  are  in  the 

same  form  as  Wedemeyer's  (28a-b) , but  their  meaning  is  different. 

These  equations  are  decoupled  by  introducing 


A = + iu^ 


and  solved  assuming  the  same  periodicity,  ei('Ct‘0'),  as  for  F ; i.e.  , m=l. 
The  results  for  and  are 


ux  = (i/2)  [(vq  + iuQ  ) e'ay  - (v  - iuQ  ) e"6y]  (A. 4a) 

00  OO  CO  oo 

v1  = - (1/2)  [(VQ  + iuQ  ) e"ay  + (VQ  - iuQ  ) e-f3y]  (A. 4b) 

OO  oo  oo  oo 
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where 


a = - Re*2  [(3  - CR)2  + eial 

a1  = (1/2)  [{2  - H (Cj)>  ir  + tan-1  {(3  - CR)/Cj}] 

6 = - Re%  [(1  + CR)2  + C:2]%  eiJ2  (A. 5) 

a2  = (1/2)  [{2  + H (Cj)>  7T  - tan’1  {(1  + CR)/Cj}] 


H CCj)  = 0 for  C < 0,  =1  for  Cj  > 0. 

The  principal  values  for  tan-1  are  used. 

Next  the  displacement  surface,  y = A(r  0,  t) , is  found.  This  is 

an  impermeable  surface  in  the  o flow;  i.e. , w = 0 on  it.  A satisfies 

a first  order  partial  differential  equation,  as  for  a steady  3-D 
boundary  layer.  For  unsteady  incompressible  flow  an  additional  term, 
rA  , appears.  The  equation  for  A is 

tr  u0oo  O - 6r)lr  + (v0oo  O ' 60)]q  + - At  = 0 (A. 6) 

where 

CO  cc 

uo»  ^ 5 / Oqoo  - u)  dy  = - / u dy 
0 0 

= (1/2)  (a'1  + g-1)  uqoo  - (i/2)  (a-1  - g"1)  vqot  , 

00  oo 

V 6°  = / Cv  - V)  dy  = - / V dy 
0 0 

= (1/2)  (a-1  + 0-1)  voco  + (i/2)  (C6-1  - g-1)  uQOT 

using  (A.4a-b). 

After  some  analysis,  employing  the  order  of  magnitude  estimates 
for  Auq  and  Avq  given  in  the  paragraph  preceding  (A.3a-d),  (A. 6)  can 

be  written  as 
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(A. 7) 


[ru^  (A  - 6)]r  + 


[v 

L o°° 


(A  - 6)]Q  + rAt 


0 


where 


26  = a-1  + B"1  - 2 (a"1  - B_1)/(l  - C)  . (A. 8) 

The  only  periodic  non-singular  solution  of  (A. 7)  is  the  constant 

A = 6 . 


From  (A. 5),  A depends  only  on  Re  and  c,  and  is  easily  found  after  the 

eigenvalue  C = C + iCT  is  determined  for  the  problem  o.  Putting  CT  = 0 
K I I 

in  (A. 5)  and  (A. 8)  gives  6 = 6c  = Wedemeyer's  displacement  thickness 
(non-dimensionalized) . Wedemeyer  also  determines  a displacement  thick- 
ness for  the  sidewall  boundary  layer.  In  our  problem  this  is  not 
necessary  since  the  proper  boundary  conditions  are  satisfied  at  the 
sidewall . 

The  displacement  thickness  is  a complex  constant.  The  interpreta- 
tion of  the  complex  A is  given  in  Reference  1.  Operationally  it  is  used 
in  the  same  manner  as  a real  displacement  thickness:  2c  is  replaced  by 
2c  - 2A  and  the  o problem  is  re-solved. 
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APPENDIX  B:  POWER  SERIES  COEFFICIENTS 


In  accordance  with  Section  III.B,  analytical  solutions  to  (9)  are 
found  for  the  interval  0 < r < e.  These  are  obtained  by  substituting 
power  series  for  u,  v,  w,  and  p into  (9)  and  determining  the  coeffi- 
cients . 


The  solution 

for  the  azimuthal 

mode  m = 

1 is  given  by 

u = E 

a.  rj  , 

v = 

E b . r^ 

o 

11 

■1 — i 

J 

j=0  3 

w = E 

d.  r^  , 

A 

P = 

E e . r^ 

o 

11 

•n 

J 

j=0  3 

The  coefficients  are  given  by  the  following  sequence  of  formulas,  where 
b , d^,  and  e^  are  arbitrary  complex  constants: 

a = i b , b=b  !,  d=0  , e =0 

o o o o o o 

a:  = 0 , bx  = 0 , d:  = d:  , ei  = ei 

H = Re-1  K2  + i C 

a2  = - (1/4)  Kd]+  (1/8)  Re  [bQ  (iH  - 1)  + e^ 

b2  = - i (1/4)  K dj_  - (3i/8)  Re  [bQ  (iH  - 1)  + . 

For  j =0,  1,  2,  ... 

Nj.  = Re  [H  aj+1  - (2  bj+1  ♦ iaj+1)] 

N2.  = Re  [H  b.+1  * (2  a.+1  - ibj+1)] 

N,  = [NJj  - i (j  + 2)  N2j]/[(j  * l)(j  + 3)] 
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d.+2  = Re  [(H  - i)  d.  - K e.]/[(j  + 1)0  + 
aj+3  a - P*j  + (3  + 4)  K dj+2]/ [(j  + 5)(j 
bj+3  = i [K  dj+2  + (j  + 4)  N.]/[(j  + 5)(j 
ej+2  = [2  aj+3  + i (j2  ^ 6j  + 7)  b.+3  - i 


3)] 

+ 3)] 

+ 3)] 
N2j]/Re 
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